!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
program data_manage
!==================================================================================================================================|

   use mpi
   implicit none

   include 'data_manage_precision.include'
   include 'data_manage_parameter.include'
!   include 'silo.inc'

!----------------------------------------------------------------------------------------------------------------------------------|
!  ro, gp, mp, te, beta, vx, vy, vz, vv, bx, by, bz, bv, jx, jy, jz, jv, fgx, fgy, fgz, fgv, fggx, fggy, fggz, fggv, 
!  fmgx, fmgy, fmgz, fmgv, fmtx, fmty, fmtz, fmtv
!----------------------------------------------------------------------------------------------------------------------------------|
   integer(4), parameter :: n_var_total = 33
   integer(4), parameter :: var_total_dim(n_var_total) = (/ 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 3, 1, 1, 1, 3, 1, 1, 1, 3,         &
                                                            1, 1, 1, 3, 1, 1, 1, 3, 1, 1, 1, 3 /)
   character(128), parameter :: var_total_name_list(n_var_total) =                                                                 &
      (/ 'ro', 'gp', 'mp', 'te', 'beta', 'vx', 'vy', 'vz', 'vv', 'bx', 'by', 'bz', 'bv', 'jx', 'jy', 'jz', 'jv',                   &
      'fgx', 'fgy', 'fgz', 'fgv', 'fggx', 'fggy', 'fggz', 'fggv', 'fmgx', 'fmgy', 'fmgz', 'fmgv', 'fmtx', 'fmty', 'fmtz', 'fmtv' /)

   character(4) :: ti_char, cpui_char
   character(128) :: datavis_path, var_name_list, directory, datavis_name
   character(256) :: header !, block_id_char
   character(128), allocatable :: var_name(:) !, mesh_name_silo(:) , var_name_silo(:, :)

   integer(4) :: i_var, j_var, n_var, len_str, format, output_method, output_level, i_level
   integer(4) :: i, j, k
   integer(4) :: i_c, j_c, k_c
   integer(4) :: i_p, j_p, k_p
   integer(4) :: i_e, j_e, k_e
   integer(4) :: ip1, jp1, kp1
   integer(4) :: nx_tmp, ny_tmp, nz_tmp
   integer(4) :: nx_len, ny_len, nz_len
   integer(4) :: nx_lv, ny_lv, nz_lv
   integer(4) :: nx_no_bnd, nx_block, nx_total
   integer(4) :: ny_no_bnd, ny_block, ny_total
   integer(4) :: nz_no_bnd, nz_block, nz_total
   integer(4) :: ini_block_x, ini_block_y, ini_block_z
   integer(4) :: ti, ci, cpui, i_block
   integer(4) :: cpu_number, ng, ratio, n_level, nc
   integer(4) :: nt, n_blocks, lv, n_dims, sys
   integer(4) :: datamap_unit, datavis_unit, datavtu_unit !, datasilo_unit
   integer(4) :: error, cpu_id, cpu_n, n_out
   integer(4) :: output_time_start, output_time_end
!   integer(4) :: optlist_id, silo_status
   integer(4), allocatable :: id(:, :)
!   integer(4), allocatable :: mesh_name_len_silo(:), mesh_type_silo(:), var_name_len_silo(:, :), var_type_silo(:, :)
!   integer(4), allocatable :: level_id_silo(:), level_data_silo(:, :)
!   integer(4), allocatable :: level_data_id_silo(:), level_type_silo(:), level_index_silo(:)
   integer(4), allocatable :: p_pos(:), n_blocks_all(:)
   integer(4), allocatable :: var_dim(:), level_len(:)
   integer(8) :: data_offset

   logical(4) :: var_flag(n_var_total)

   real(8) :: dx, dy, dz
   real(8) :: rangex, rangey, rangez
   real(8) :: zero_posx, zero_posy, zero_posz
   real(8) :: gm, cfl, gmm1, eps
   real(8), dimension(:), allocatable :: t, dt
   real(8), dimension(:, :), allocatable :: x_arr, y_arr, z_arr
   real(8), dimension(:), allocatable :: x, y, z, xh, yh, zh, x_lv, y_lv, z_lv
!----------------------------------------------------------------------------------------------------------------------------------|
!  ro, gp, mp, te, beta, vx, vy, vz, vv, bx, by, bz, bv, jx, jy, jz, jv, fgx, fgy, fgz, fgv, fggx, fggy, fggz, fggv, 
!  fmgx, fmgy, fmgz, fmgv, fmtx, fmty, fmtz, fmtv
!----------------------------------------------------------------------------------------------------------------------------------|
   real(output_real8_flag), dimension(:, :, :), allocatable :: ro, mx, my, mz, bx, by, bz, en, gp, gx, gy, gz,                     &
                                                               roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh,            &
                                                               jxh, jyh, jzh, fgxh, fgyh, fgzh, fggxh, fggyh, fggzh,               &
                                                               fmgxh, fmgyh, fmgzh, fmtxh, fmtyh, fmtzh,                           &
                                                               ro_lv,vx_lv, vy_lv, vz_lv, bx_lv, by_lv, bz_lv, gp_lv

   real(output_real8_flag), allocatable :: var(:, :, :, :, :), var_max(:), var_min(:)
   namelist /datavis_options/                                                                                                      &
      output_method,                                                                                                               &
      output_level,                                                                                                                &
      output_time_start,                                                                                                           &
      output_time_end,                                                                                                             &
      datavis_path,                                                                                                                &
      var_name_list,                                                                                                               &
      format,                                                                                                                      &
      header

!----------------------------------------------------------------------------------------------------------------------------------|
!  start
!----------------------------------------------------------------------------------------------------------------------------------|
   call mpi_init(error)
   call mpi_comm_size(mpi_comm_world, cpu_n, error)
   call mpi_comm_rank(mpi_comm_world, cpu_id, error)

   if (cpu_id .eq. 0) then
      write(*, *) 
      write(*, '(a)') '|------------------|'
      write(*, '(a)') ' Data manage starts:'
      write(*, '(a)') '|------------------|'
      write(*, *)
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  check no reduplicate variable name
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, n_var_total
      i_var = 0
      do i = 1, n_var_total
         if (index(trim(var_total_name_list(j)), trim(var_total_name_list(i))) .ne. 0) i_var = i_var + 1
      enddo
      if (i_var .gt. 1) then
         write(*, '(a, 1x, a)') ' ERROR: reduplicate variable name:', var_total_name_list(j)
         do i = 1, n_var_total
            if (index(trim(var_total_name_list(j)), trim(var_total_name_list(i))) .ne. 0)                                          &
               write(*, '(1x, a, 1x, a, 1x, a)'), trim(var_total_name_list(j)), '<-', trim(var_total_name_list(i))
         enddo
         write(*, *)
         stop
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  prepare
!----------------------------------------------------------------------------------------------------------------------------------|
   eps = 1.0d-12
   datamap_unit = 10
   datavis_unit = 11
   datavtu_unit = 12
!   datasilo_unit = 13

   open (datamap_unit, file = trim(parameters_filename))
   read(datamap_unit, nml = datavis_options)
   close (datamap_unit)

   call system("mkdir -p " // trim(datavis_path))

   var_flag = .false.

   len_str = len(var_name_list)
   do i = 1, len_str
      if (var_name_list(i:i) .ge. 'A' .and. var_name_list(i:i) .le. 'Z')                                                           &
         var_name_list(i:i) = char(ichar(var_name_list(i:i)) - ichar('A') +  ichar('a'))
   enddo

   n_var = 0
   do i_var = 1, n_var_total
      if (index(var_name_list, trim(var_total_name_list(i_var))) .ne. 0) then
         var_flag(i_var) = .true.
         n_var = n_var + 1
      endif
   enddo

   if (cpu_id .eq. 0) then
      write(*, '(a, i4)') ' Total number of variables is:', n_var
      write(*, '(a)') ' The following variables will be outputted for visualization:'
      write(*, *)
   endif

   allocate (var_name(n_var), var_dim(n_var), var_max(n_var), var_min(n_var))

   if (cpu_id .eq. 0) then
      do i_var = 1, n_var_total
         if (var_flag(i_var)) write(*, '(1x, a8, 1x)', advance = 'no') trim(var_total_name_list(i_var))
      enddo
      write(*, *)
   endif

   var_name = ''
   var_dim = 0

   j_var = 1
   do i_var = 1, n_var_total
      if (var_flag(i_var)) then
         var_name(j_var) = var_total_name_list(i_var)
         var_dim(j_var) = var_total_dim(i_var)
         j_var = j_var + 1
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  read parameters
!----------------------------------------------------------------------------------------------------------------------------------|
   open(datamap_unit, file = trim(datamap_path) // 'params.mapout')
   read(datamap_unit, '(a)') directory
   read(datamap_unit, '(t32, i)') cpu_number
   read(datamap_unit, '(t32, i)') n_dims
   read(datamap_unit, '(t32, i)') sys
   read(datamap_unit, '(t32, i)') ng
   read(datamap_unit, '(t32, i)') ini_block_x
   read(datamap_unit, '(t32, i)') ini_block_y
   read(datamap_unit, '(t32, i)') ini_block_z
   read(datamap_unit, '(t32, i)') nx_no_bnd
   read(datamap_unit, '(t32, i)') nx_block
   read(datamap_unit, '(t32, i)') nx_total
   read(datamap_unit, '(t32, i)') ny_no_bnd
   read(datamap_unit, '(t32, i)') ny_block
   read(datamap_unit, '(t32, i)') ny_total
   read(datamap_unit, '(t32, i)') nz_no_bnd
   read(datamap_unit, '(t32, i)') nz_block
   read(datamap_unit, '(t32, i)') nz_total
   read(datamap_unit, '(t32, i)') ratio
   read(datamap_unit, '(t32, i)') n_level
   read(datamap_unit, '(t32, es)') gm
   read(datamap_unit, '(t32, es)') rangex
   read(datamap_unit, '(t32, es)') rangey
   read(datamap_unit, '(t32, es)') rangez
   read(datamap_unit, '(t32, es)') zero_posx
   read(datamap_unit, '(t32, es)') zero_posy
   read(datamap_unit, '(t32, es)') zero_posz
   read(datamap_unit, '(t32, es)') cfl

   if (cpu_id .eq. 0) then
      write(*, *)
      write(*, '(a, 1x, a)') ' Work space:', trim(directory)
      write(*, *)
      write(*, "(a, i5)") ' Base x grid points=', (nx_total - ng * 2)
      write(*, "(a, i5)") ' Base y grid points=', (ny_total - ng * 2)
      write(*, "(a, i5)") ' Base z grid points=', (nz_total - ng * 2)
      write(*, *)
      write(*, '(a, i)') ' Output method:', output_method
      write(*, *)
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  read the number of times of outputs
!----------------------------------------------------------------------------------------------------------------------------------|
   open (datamap_unit, file = trim(datamap_path) // 'nt.mapout', form = 'unformatted')
   read(datamap_unit) nt
   close (datamap_unit)

!----------------------------------------------------------------------------------------------------------------------------------|
!  ro, gp, mp, te, beta, vx, vy, vz, vv, bx, by, bz, bv, jx, jy, jz, jv, fgx, fgy, fgz, fgv, fggx, fggy, fggz, fggv, 
!  fmgx, fmgy, fmgz, fmgv, fmtx, fmty, fmtz, fmtv
!----------------------------------------------------------------------------------------------------------------------------------|
   nx_lv = (nx_total - ng * 2) * 2 ** (output_level - 1) + 2 * ng
   ny_lv = (ny_total - ng * 2) * 2 ** (output_level - 1) + 2 * ng
   nz_lv = (nz_total - ng * 2) * 2 ** (output_level - 1) + 2 * ng
   allocate(t(nt), dt(nt), n_blocks_all(nt), p_pos(n_dims), id(n_dims, n_level), level_len(n_level),                               &
      x(nx_block), y(ny_block), z(nz_block),                                                                                       &
      xh(nx_block), yh(ny_block), zh(nz_block),                                                                                    &
      x_lv(nx_lv), y_lv(ny_lv), z_lv(nz_lv),                                                                                       &
      ro(nx_block, ny_block, nz_block), mx(nx_block, ny_block, nz_block), my(nx_block, ny_block, nz_block),                        &
      mz(nx_block, ny_block, nz_block), bx(nx_block, ny_block, nz_block), by(nx_block, ny_block, nz_block),                        &
      bz(nx_block, ny_block, nz_block), en(nx_block, ny_block, nz_block), gp(nx_block, ny_block, nz_block),                        &
      gx(nx_block, ny_block, nz_block), gy(nx_block, ny_block, nz_block), gz(nx_block, ny_block, nz_block),                        &
      roh(nx_block, ny_block, nz_block), gph(nx_block, ny_block, nz_block), mph(nx_block, ny_block, nz_block),                     &
      teh(nx_block, ny_block, nz_block), betah(nx_block, ny_block, nz_block), vxh(nx_block, ny_block, nz_block),                   &
      vyh(nx_block, ny_block, nz_block), vzh(nx_block, ny_block, nz_block), bxh(nx_block, ny_block, nz_block),                     &
      byh(nx_block, ny_block, nz_block), bzh(nx_block, ny_block, nz_block), jxh(nx_block, ny_block, nz_block),                     &
      jyh(nx_block, ny_block, nz_block), jzh(nx_block, ny_block, nz_block), fgxh(nx_block, ny_block, nz_block),                    &
      fgyh(nx_block, ny_block, nz_block), fgzh(nx_block, ny_block, nz_block), fggxh(nx_block, ny_block, nz_block),                 &
      fggyh(nx_block, ny_block, nz_block), fggzh(nx_block, ny_block, nz_block), fmgxh(nx_block, ny_block, nz_block),               &
      fmgyh(nx_block, ny_block, nz_block), fmgzh(nx_block, ny_block, nz_block), fmtxh(nx_block, ny_block, nz_block),               &
      fmtyh(nx_block, ny_block, nz_block), fmtzh(nx_block, ny_block, nz_block),                                                    &
      ro_lv(nx_lv, ny_lv, nz_lv), gp_lv(nx_lv, ny_lv, nz_lv),                                                                      &
      vx_lv(nx_lv, ny_lv, nz_lv), vy_lv(nx_lv, ny_lv, nz_lv), vz_lv(nx_lv, ny_lv, nz_lv),                                          &
      bx_lv(nx_lv, ny_lv, nz_lv), by_lv(nx_lv, ny_lv, nz_lv), bz_lv(nx_lv, ny_lv, nz_lv))

!----------------------------------------------------------------------------------------------------------------------------------|
   gmm1 = gm - 1.0d0

   n_blocks_all = 0

   output_time_start = max(0, output_time_start) + 1
   output_time_end = min(nt - 1, output_time_end) + 1

   n_out = 0

   call mpi_barrier(mpi_comm_world, error)

   do ti = output_time_start, output_time_end
!   do ti = 1, nt

      if (ti .ne. output_time_start + cpu_id + cpu_n * n_out) cycle
      n_out = n_out + 1

      write(ti_char, '(i4.4)') ti - 1

!----------------------------------------------------------------------------------------------------------------------------------|
!     read t, dt
!----------------------------------------------------------------------------------------------------------------------------------|
      open(datamap_unit, file = trim(datamap_path) // 'info.t' // ti_char // '.mapout', form = 'unformatted')
      read(datamap_unit) t(ti)
      read(datamap_unit) dt(ti)
      close(datamap_unit)

      datavis_name = trim(datavis_path) // 'data' // ti_char

      write(*, '(a, es10.3)') ' Time =', t(ti)

      var_min = 1.0e10
      var_max = -1.0e10
      i_block = 0

      if (output_method .eq. 2) open(datavis_unit, file = trim(datavis_name) // '.mapvis', form = 'unformatted')

      if (output_method .eq. 3) then
         data_offset = 0
         if (format .eq. 0) then
            open(datavtu_unit, file = trim(datavis_name) // '.vtu', form = 'binary')
            write(datavtu_unit) '<?xml version="1.0"?>' // char(10)
            write(datavtu_unit) '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">' // char(10)
            write(datavtu_unit) '<UnstructuredGrid>' // char(10)
         else
            open(datavtu_unit, file = trim(datavis_name) // '.vtu')
            write(datavtu_unit, '(a)') '<?xml version="1.0"?>'
            write(datavtu_unit, '(a)') '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">'
            write(datavtu_unit, '(a)') '<UnstructuredGrid>'
         endif

         if (format .eq. 0) then
            do cpui = 0, cpu_number - 1
               write(cpui_char, '(i4.4)') cpui

               open(datamap_unit, file = trim(datamap_path) // 'data.t' // ti_char // '.cpu' // cpui_char // '.mapout',            &
                  form = 'unformatted')
               read(datamap_unit) n_blocks_all(ti)
               read(datamap_unit) level_len
               read(datamap_unit) n_blocks
               do ci = 1, n_blocks
                  read(datamap_unit) id
                  read(datamap_unit) lv, nc, p_pos
                  read(datamap_unit) dx
                  read(datamap_unit) dy
                  read(datamap_unit) dz

                  read(datamap_unit) x
                  read(datamap_unit) y
                  read(datamap_unit) z
                  read(datamap_unit) ro
                  read(datamap_unit) mx
                  read(datamap_unit) my
                  read(datamap_unit) mz
                  read(datamap_unit) bx
                  read(datamap_unit) by
                  read(datamap_unit) bz
                  read(datamap_unit) en
                  read(datamap_unit) gx
                  read(datamap_unit) gy
                  read(datamap_unit) gz
                  if (nc .ne. 0) cycle
                  call write_vtu_header(datavtu_unit, var_name,                                                                    &
                     nx_block - 2 * ng + 1, ny_block - 2 * ng + 1, nz_block - 2 * ng + 1, n_var, var_dim, format, data_offset)
               enddo
               close(datamap_unit)
            enddo
         endif
      endif

      if (output_method .eq. 3) then
         if (format .eq. 0) then
            write(datavtu_unit) '</UnstructuredGrid>' // char(10)
            write(datavtu_unit) '<AppendedData encoding="raw">' // char(10)
            write(datavtu_unit) '_'
         endif
      endif

!      if (output_method .eq. 4) then
!         error = dbcreate(trim(datavis_name) // '.silo', len(trim(datavis_name) // '.silo'),                                       &
!            db_clobber, db_local, "AMR data", 8, db_pdb, datasilo_unit)
!      endif

      do cpui = 0, cpu_number - 1
!----------------------------------------------------------------------------------------------------------------------------------|
!        read n_blocks
!----------------------------------------------------------------------------------------------------------------------------------|
         write(cpui_char, '(i4.4)') cpui

         open(datamap_unit, file = trim(datamap_path) // 'data.t' // ti_char // '.cpu' // cpui_char // '.mapout',                  &
            form = 'unformatted')
         read(datamap_unit) n_blocks_all(ti)
         read(datamap_unit) level_len
         read(datamap_unit) n_blocks

         if (cpui .eq. 0) then
            if (output_method .eq. 1) then
               allocate(var(sum(var_dim), nx_block, ny_block, nz_block, n_blocks_all(ti)),                                         &
                                      x_arr(nx_block, n_blocks_all(ti)),                                                           &
                                      y_arr(ny_block, n_blocks_all(ti)),                                                           &
                                      z_arr(nz_block, n_blocks_all(ti)))
            endif
            if (output_method .eq. 3) then
               allocate(var(sum(var_dim), nx_block, ny_block, nz_block, 1),                                                        &
                                      x_arr(nx_block, 1),                                                                          &
                                      y_arr(ny_block, 1),                                                                          &
                                      z_arr(nz_block, 1))
            endif
!            if (output_method .eq. 4) then
!               allocate(var(sum(var_dim), nx_block, ny_block, nz_block, 1),                                                        &
!                                      x_arr(nx_block, 1),                                                                          &
!                                      y_arr(ny_block, 1),                                                                          &
!                                      z_arr(nz_block, 1),                                                                          &
!                                      mesh_name_silo(n_blocks_all(ti)),                                                            &
!                                      mesh_name_len_silo(n_blocks_all(ti)),                                                        &
!                                      mesh_type_silo(n_blocks_all(ti)),                                                            &
!                                      var_name_silo(n_var, n_blocks_all(ti)),                                                      &
!                                      var_name_len_silo(n_var, n_blocks_all(ti)),                                                  &
!                                      var_type_silo(n_var, n_blocks_all(ti)),                                                      &
!                                      level_id_silo(n_level),                                                                      &
!                                      level_data_silo(n_level, maxval(level_len)),                                                 &
!                                      level_type_silo(n_level),                                                                    &
!                                      level_index_silo(n_level),                                                                   &
!                                      level_data_id_silo(n_level))
!               level_index_silo = 0
!            endif
         endif

         do ci = 1, n_blocks
!----------------------------------------------------------------------------------------------------------------------------------|
!           read x, ro, mx, my, mz, bx, by, bz, en, gx, gy, gz
!----------------------------------------------------------------------------------------------------------------------------------|
            read(datamap_unit) id
            read(datamap_unit) lv, nc, p_pos
            read(datamap_unit) dx
            read(datamap_unit) dy
            read(datamap_unit) dz

            read(datamap_unit) x
            read(datamap_unit) y
            read(datamap_unit) z
            read(datamap_unit) ro
            read(datamap_unit) mx
            read(datamap_unit) my
            read(datamap_unit) mz
            read(datamap_unit) bx
            read(datamap_unit) by
            read(datamap_unit) bz
            read(datamap_unit) en
            read(datamap_unit) gx
            read(datamap_unit) gy
            read(datamap_unit) gz

            i_e = nx_block - 1
            j_e = ny_block - 1
            k_e = nz_block - 1

            gp = (en - (mx ** 2 + my ** 2 + mz ** 2) / ro / 2.0 - (bx ** 2 + by ** 2 + bz ** 2) / 2.0) * gmm1

            do i = 1, i_e
               xh(i) = (x(i) + x(i + 1)) / 2.0d0
            enddo

            do j = 1, j_e
               yh(j) = (y(j) + y(j + 1)) / 2.0d0
            enddo

            do k = 1, k_e
               zh(k) = (z(k) + z(k + 1)) / 2.0d0
            enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!           ro, gp, mp, te, beta, vx, vy, vz, vv, bx, by, bz, bv, jx, jy, jz, jv, fgx, fgy, fgz, fgv, fggx, fggy, fggz, fggv, 
!           fmgx, fmgy, fmgz, fmgv, fmtx, fmty, fmtz, fmtv
!----------------------------------------------------------------------------------------------------------------------------------|
            do k = 1, k_e
               kp1 = k + 1
               do j = 1, j_e
                  jp1 = j + 1
                  do i = 1, i_e
                     ip1 = i + 1
                     roh(i, j, k) = (ro(i, j, k) + ro(ip1, j, k) + ro(i, jp1, k) + ro(ip1, jp1, k) +                               &
                                  ro(i, j, kp1) + ro(ip1, j, kp1) + ro(i, jp1, kp1) + ro(ip1, jp1, kp1)) / 8.0
                     vxh(i, j, k) = (mx(i, j, k) + mx(ip1, j, k) + mx(i, jp1, k) + mx(ip1, jp1, k) +                               &
                                  mx(i, j, kp1) + mx(ip1, j, kp1) + mx(i, jp1, kp1) + mx(ip1, jp1, kp1)) / 8.0 / roh(i, j, k)
                     vyh(i, j, k) = (my(i, j, k) + my(ip1, j, k) + my(i, jp1, k) + my(ip1, jp1, k) +                               &
                                  my(i, j, kp1) + my(ip1, j, kp1) + my(i, jp1, kp1) + my(ip1, jp1, kp1)) / 8.0 / roh(i, j, k)
                     vzh(i, j, k) = (mz(i, j, k) + mz(ip1, j, k) + mz(i, jp1, k) + mz(ip1, jp1, k) +                               &
                                  mz(i, j, kp1) + mz(ip1, j, kp1) + mz(i, jp1, kp1) + mz(ip1, jp1, kp1)) / 8.0 / roh(i, j, k)
                     bxh(i, j, k) = (bx(i, j, k) + bx(ip1, j, k) + bx(i, jp1, k) + bx(ip1, jp1, k) +                               &
                                  bx(i, j, kp1) + bx(ip1, j, kp1) + bx(i, jp1, kp1) + bx(ip1, jp1, kp1)) / 8.0
                     byh(i, j, k) = (by(i, j, k) + by(ip1, j, k) + by(i, jp1, k) + by(ip1, jp1, k) +                               &
                                  by(i, j, kp1) + by(ip1, j, kp1) + by(i, jp1, kp1) + by(ip1, jp1, kp1)) / 8.0
                     bzh(i, j, k) = (bz(i, j, k) + bz(ip1, j, k) + bz(i, jp1, k) + bz(ip1, jp1, k) +                               &
                                  bz(i, j, kp1) + bz(ip1, j, kp1) + bz(i, jp1, kp1) + bz(ip1, jp1, kp1)) / 8.0
                     gph(i, j, k) = (gp(i, j, k) + gp(ip1, j, k) + gp(i, jp1, k) + gp(ip1, jp1, k) +                               &
                                  gp(i, j, kp1) + gp(ip1, j, kp1) + gp(i, jp1, kp1) + gp(ip1, jp1, kp1)) / 8.0
                     mph(i, j, k) = (bxh(i, j, k) * bxh(i, j, k) + byh(i, j, k) * byh(i, j, k) + bzh(i, j, k) * bzh(i, j, k)) / 2.0
                     teh(i, j, k) = gph(i, j, k) / roh(i, j, k)
                     betah(i, j, k) = gph(i, j, k) / mph(i, j, k)
                     fgxh(i, j, k) = (gx(i, j, k) + gx(ip1, j, k) + gx(i, jp1, k) + gx(ip1, jp1, k) +                              &
                                  gx(i, j, kp1) + gx(ip1, j, kp1) + gx(i, jp1, kp1) + gx(ip1, jp1, kp1)) / 8.0 * roh(i, j, k)
                     fgyh(i, j, k) = (gy(i, j, k) + gy(ip1, j, k) + gy(i, jp1, k) + gy(ip1, jp1, k) +                              &
                                  gy(i, j, kp1) + gy(ip1, j, kp1) + gy(i, jp1, kp1) + gy(ip1, jp1, kp1)) / 8.0 * roh(i, j, k)
                     fgzh(i, j, k) = (gz(i, j, k) + gz(ip1, j, k) + gz(i, jp1, k) + gz(ip1, jp1, k) +                              &
                                  gz(i, j, kp1) + gz(ip1, j, kp1) + gz(i, jp1, kp1) + gz(ip1, jp1, kp1)) / 8.0 * roh(i, j, k)
                  enddo
               enddo
            enddo

            do k = 2, k_e - 1
               do j = 2, j_e - 1
                  do i = 2, i_e - 1
                     jxh(i, j, k) = (bzh(i, j + 1, k) - bzh(i, j - 1, k)) / dy / 2.0 -                                             &
                                    (byh(i, j, k + 1) - byh(i, j, k - 1)) / dz / 2.0
                     jyh(i, j, k) = (bxh(i, j, k + 1) - bxh(i, j, k - 1)) / dz / 2.0 -                                             &
                                    (bzh(i + 1, j, k) - bzh(i - 1, j, k)) / dx / 2.0
                     jzh(i, j, k) = (byh(i + 1, j, k) - byh(i - 1, j, k)) / dx / 2.0 -                                             &
                                    (bxh(i, j + 1, k) - bxh(i, j - 1, k)) / dy / 2.0
                     fggxh(i, j, k) = (gph(i + 1, j, k) - gph(i - 1, j, k)) / dx / 2.0
                     fggyh(i, j, k) = (gph(i, j + 1, k) - gph(i, j - 1, k)) / dy / 2.0
                     fggzh(i, j, k) = (gph(i, j, k + 1) - gph(i, j, k - 1)) / dz / 2.0
                     fmgxh(i, j, k) = (mph(i + 1, j, k) - mph(i - 1, j, k)) / dx / 2.0
                     fmgyh(i, j, k) = (mph(i, j + 1, k) - mph(i, j - 1, k)) / dy / 2.0
                     fmgzh(i, j, k) = (mph(i, j, k + 1) - mph(i, j, k - 1)) / dz / 2.0
                     fmtxh(i, j, k) = bxh(i, j, k) * (bxh(i + 1, j, k) - bxh(i - 1, j, k)) / dx / 2.0 +                            &
                                      byh(i, j, k) * (bxh(i, j + 1, k) - bxh(i, j - 1, k)) / dy / 2.0 +                            &
                                      bzh(i, j, k) * (bxh(i, j, k + 1) - bxh(i, j, k - 1)) / dz / 2.0
                     fmtyh(i, j, k) = bxh(i, j, k) * (byh(i + 1, j, k) - byh(i - 1, j, k)) / dx / 2.0 +                            &
                                      byh(i, j, k) * (byh(i, j + 1, k) - byh(i, j - 1, k)) / dy / 2.0 +                            &
                                      bzh(i, j, k) * (byh(i, j, k + 1) - byh(i, j, k - 1)) / dz / 2.0
                     fmtzh(i, j, k) = bxh(i, j, k) * (bzh(i + 1, j, k) - bzh(i - 1, j, k)) / dx / 2.0 +                            &
                                      byh(i, j, k) * (bzh(i, j + 1, k) - bzh(i, j - 1, k)) / dy / 2.0 +                            &
                                      bzh(i, j, k) * (bzh(i, j, k + 1) - bzh(i, j, k - 1)) / dz / 2.0
                  enddo
               enddo
            enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!           idl base level
!----------------------------------------------------------------------------------------------------------------------------------|
            if (output_method .eq. 2) then
               if (lv .eq. output_level) then
                  nx_tmp = 0
                  ny_tmp = 0
                  nz_tmp = 0
                  do i_level = 1, output_level
                     nx_tmp = nx_tmp + (id(1, i_level) - 1) * nx_no_bnd * 2 ** (output_level - i_level)
                     ny_tmp = ny_tmp + (id(2, i_level) - 1) * ny_no_bnd * 2 ** (output_level - i_level)
                     nz_tmp = nz_tmp + (id(3, i_level) - 1) * nz_no_bnd * 2 ** (output_level - i_level)
                  enddo
                  nx_tmp = nx_tmp + ng
                  ny_tmp = ny_tmp + ng
                  nz_tmp = nz_tmp + ng

                  do i = 0, nx_no_bnd
                     i_c = i + ng
                     i_p = i + nx_tmp
                     x_lv(i_p) = xh(i_c)
                  enddo

                  do j = 0, ny_no_bnd
                     j_c = j + ng
                     j_p = j + ny_tmp
                     y_lv(j_p) = yh(j_c)
                  enddo

                  do k = 0, nz_no_bnd
                     k_c = k + ng
                     k_p = k + nz_tmp
                     z_lv(k_p) = zh(k_c)
                  enddo

                  do k = 0, nz_no_bnd
                     k_c = k + ng
                     k_p = k + nz_tmp
                     do j = 0, ny_no_bnd
                        j_c = j + ng
                        j_p = j + ny_tmp
                        do i = 0, nx_no_bnd
                           i_c = i + ng
                           i_p = i + nx_tmp
                           ro_lv(i_p, j_p, k_p) = roh(i_c, j_c, k_c)
                           vx_lv(i_p, j_p, k_p) = vxh(i_c, j_c, k_c)
                           vy_lv(i_p, j_p, k_p) = vyh(i_c, j_c, k_c)
                           vz_lv(i_p, j_p, k_p) = vzh(i_c, j_c, k_c)
                           bx_lv(i_p, j_p, k_p) = bxh(i_c, j_c, k_c)
                           by_lv(i_p, j_p, k_p) = byh(i_c, j_c, k_c)
                           bz_lv(i_p, j_p, k_p) = bzh(i_c, j_c, k_c)
                           gp_lv(i_p, j_p, k_p) = gph(i_c, j_c, k_c)
                        enddo
                     enddo
                  enddo
               else if (lv .lt. output_level) then
                  nx_len = nx_no_bnd * 2 ** (output_level - lv)
                  ny_len = ny_no_bnd * 2 ** (output_level - lv)
                  nz_len = nz_no_bnd * 2 ** (output_level - lv)

                  nx_tmp = 0
                  ny_tmp = 0
                  nz_tmp = 0
                  do i_level = 1, lv
                     nx_tmp = nx_tmp + (id(1, i_level) - 1) * nx_no_bnd * 2 ** (output_level - i_level)
                     ny_tmp = ny_tmp + (id(2, i_level) - 1) * ny_no_bnd * 2 ** (output_level - i_level)
                     nz_tmp = nz_tmp + (id(3, i_level) - 1) * nz_no_bnd * 2 ** (output_level - i_level)
                  enddo

!                  do i = 0, nx_len
!                     i_c = i + ng
!                     i_p = i + nx_tmp
!                     if () then
!                     x_lv(i_p) = xh(i_c)
!                  enddo

               endif
            endif

            if (nc .ne. 0) cycle

!----------------------------------------------------------------------------------------------------------------------------------|
!           calculate the maximum and minimum value
!----------------------------------------------------------------------------------------------------------------------------------|
            call extremum_calculate(var_min, var_max, roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh, jxh, jyh, jzh,      &
               fgxh, fgyh, fgzh, fggxh, fggyh, fggzh, fmgxh, fmgyh, fmgzh, fmtxh, fmtyh, fmtzh, var_flag,                          &
               n_var, nx_block, ny_block, nz_block, ng, n_var_total)

!----------------------------------------------------------------------------------------------------------------------------------|
!           var array and coordinate array
!----------------------------------------------------------------------------------------------------------------------------------|
            if (output_method .eq. 1) then
               i_block = i_block + 1
               call data_collect(var, roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh, jxh, jyh, jzh,                      &
                  fgxh, fgyh, fgzh, fggxh, fggyh, fggzh, fmgxh, fmgyh, fmgzh, fmtxh, fmtyh, fmtzh, var_flag,                       &
                  x_arr, y_arr, z_arr, xh, yh, zh,                                                                                 &
                  sum(var_dim), nx_block, ny_block, nz_block, n_blocks_all(ti), n_var_total, i_block)
            endif
            if (output_method .eq. 3) then
               i_block = 1
               call data_collect(var, roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh, jxh, jyh, jzh,                      &
                  fgxh, fgyh, fgzh, fggxh, fggyh, fggzh, fmgxh, fmgyh, fmgzh, fmtxh, fmtyh, fmtzh, var_flag,                       &
                  x_arr, y_arr, z_arr, xh, yh, zh,                                                                                 &
                  sum(var_dim), nx_block, ny_block, nz_block, i_block, n_var_total, i_block)
            endif
            if (output_method .eq. 4) then
               i_block = i_block + 1
               call data_collect(var, roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh, jxh, jyh, jzh,                      &
                  fgxh, fgyh, fgzh, fggxh, fggyh, fggzh, fmgxh, fmgyh, fmgzh, fmtxh, fmtyh, fmtzh, var_flag,                       &
                  x_arr, y_arr, z_arr, xh, yh, zh,                                                                                 &
                  sum(var_dim), nx_block, ny_block, nz_block, 1, n_var_total, 1)
            endif

            if (output_method .eq. 3) then
               call write_vtu_data(datavtu_unit, var_name,                                                                         &
                  var(:, ng : nx_block - ng, ng : ny_block - ng, ng : nz_block - ng, i_block),                                     &
                  real(x_arr(ng : nx_block - ng, i_block)),                                                                        &
                  real(y_arr(ng : ny_block - ng, i_block)),                                                                        &
                  real(z_arr(ng : nz_block - ng, i_block)),                                                                        &
                  nx_block - 2 * ng + 1, ny_block - 2 * ng + 1, nz_block - 2 * ng + 1, n_var, var_dim, format)
            endif

!            if (output_method .eq. 4) then
!               write(block_id_char, '(i10.10)') i_block
!               block_id_char = '/block' // trim(block_id_char)
!               mesh_name_silo(i_block) = trim(block_id_char) // '/' // 'amr_mesh'
!               mesh_name_len_silo(i_block) = len_trim(mesh_name_silo(i_block))

!               mesh_type_silo(i_block) = db_quadmesh
!               do i_var = 1, n_var
!                  var_name_silo(i_var, i_block) = trim(block_id_char) // '/' // trim(var_name(i_var))
!                  var_name_len_silo(i_var, i_block) = len_trim(var_name_silo(i_var, i_block))
!                  var_type_silo(i_var, i_block) = db_quadvar
!               enddo

!               level_index_silo(lv) = level_index_silo(lv) + 1
!               level_data_silo(lv, level_index_silo(lv)) = i_block
!               call write_silo(datasilo_unit, var_name,                                                                            &
!                  var(:, ng : nx_block - ng, ng : ny_block - ng, ng : nz_block - ng, 1),                                           &
!                  real(x_arr(ng : nx_block - ng, 1)), real(y_arr(ng : ny_block - ng, 1)), real(z_arr(ng : nz_block - ng, 1)),      &
!                  nx_block - 2 * ng + 1, ny_block - 2 * ng + 1, nz_block - 2 * ng + 1, n_var, var_dim, block_id_char)
!            endif

!----------------------------------------------------------------------------------------------------------------------------------|
!           for the base blocks
!----------------------------------------------------------------------------------------------------------------------------------|
!            if (lv .eq. 1 .and. output_method .eq. 1) then
!               call data_merge (var, x_base, y_base, z_base,                                                                       &
!                  roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh, jxh, jyh, jzh, fgxh, fgyh, fgzh,                        &
!                  fggxh, fggyh, fggzh, fmgxh, fmgyh, fmgzh, fmtxh, fmtyh, fmtzh, xh, yh, zh, var_flag, p_pos,                      &
!                  sum(var_dim), nx_total, ny_total, nz_total, nx_block, ny_block, nz_block,                                        &
!                  nx_no_bnd, ny_no_bnd, nz_no_bnd, ng, n_var_total, n_dims)
!            endif

!            if (lv .eq. 1 .and. output_method .eq. 3) then
!               call data_merge (var, x_base, y_base, z_base,                                                                       &
!                  roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh, jxh, jyh, jzh, fgxh, fgyh, fgzh,                        &
!                  fggxh, fggyh, fggzh, fmgxh, fmgyh, fmgzh, fmtxh, fmtyh, fmtzh, xh, yh, zh, var_flag, p_pos,                      &
!                  sum(var_dim), nx_total, ny_total, nz_total, nx_block, ny_block, nz_block,                                        &
!                  nx_no_bnd, ny_no_bnd, nz_no_bnd, ng, n_var_total, n_dims)
!            endif
         enddo
         close(datamap_unit)
      enddo

!      if (output_method .eq. 4) then
!         error = dbmkoptlist(10, optlist_id)
!         error = dbaddiopt(optlist_id, dbopt_cycle, ti)
!         error = dbadddopt(optlist_id, dbopt_dtime, t(ti))
!         error = dbaddcopt(optlist_id, dbopt_mrgtree_name, 'mrgtree', 7)

!         error = dbset2dstrlen(128)
!         error = dbputmmesh(datasilo_unit, 'amr_mesh', 8, i_block,                                                                 &
!            mesh_name_silo, mesh_name_len_silo, mesh_type_silo, optlist_id, silo_status)

!         error = dbset2dstrlen(128)
!         do i_var = 1, n_var
!            error = dbputmvar(datasilo_unit, trim(var_name(i_var)), len(trim(var_name(i_var))),                                    &
!               i_block, var_name_silo(i_var, :), var_name_len_silo(i_var, :), var_type_silo(i_var, :),                             &
!               optlist_id, silo_status)
!         enddo
!         do i_level = 1, n_level
!            level_id_silo(i_level) = i_level
!            level_type_silo(i_level) = db_blockcent
!            level_data_id_silo(i_level) = dbmkptr(level_data_silo(i_level, :))
!         enddo
!!print*, mesh_name_len_silo(2)
!!         error = dbputgrplmap(datasilo_unit, "amr_levels", 10, n_level,                                                            &
!!            level_type_silo, level_len, 0, level_data_id_silo, 0, 0, 0, silo_status)
!         error = dbfreeoptlist(optlist_id)
!         error = dbset2dstrlen(32)
!         error = dbclose(datasilo_unit)
!      endif

      if (output_method .eq. 3) then
         if (format .eq. 0) then
            write(datavtu_unit) '</AppendedData>' // char(10)
            write(datavtu_unit) '</VTKFile>'
            close(datavtu_unit)
         else
            write(datavtu_unit, '(a)') '</UnstructuredGrid>'
            write(datavtu_unit, '(a)') '</VTKFile>'
            close(datavtu_unit)
         endif
      endif
      if (output_method .eq. 2) then
         write(datavis_unit) output_real8_flag
         write(datavis_unit) nx_lv - 2 * ng + 1, ny_lv - 2 * ng + 1, nz_lv - 2 * ng + 1
         write(datavis_unit) x_lv(ng : nx_lv - ng), y_lv(ng : ny_lv - ng), z_lv(ng : nz_lv - ng)
         write(datavis_unit) ro_lv(ng : nx_lv - ng, ng : ny_lv - ng, ng : nz_lv - ng)
         write(datavis_unit) vx_lv(ng : nx_lv - ng, ng : ny_lv - ng, ng : nz_lv - ng)
         write(datavis_unit) vy_lv(ng : nx_lv - ng, ng : ny_lv - ng, ng : nz_lv - ng)
         write(datavis_unit) vz_lv(ng : nx_lv - ng, ng : ny_lv - ng, ng : nz_lv - ng)
         write(datavis_unit) bx_lv(ng : nx_lv - ng, ng : ny_lv - ng, ng : nz_lv - ng)
         write(datavis_unit) by_lv(ng : nx_lv - ng, ng : ny_lv - ng, ng : nz_lv - ng)
         write(datavis_unit) bz_lv(ng : nx_lv - ng, ng : ny_lv - ng, ng : nz_lv - ng)
         write(datavis_unit) gp_lv(ng : nx_lv - ng, ng : ny_lv - ng, ng : nz_lv - ng)
         close(datavis_unit)
      endif

      if (output_method .eq. 1) then
         call write_vtk_unstructured(trim(datavis_name) // '.vtk', var_name,                                                       &
            var(:, ng : nx_block - ng, ng : ny_block - ng, ng : nz_block - ng, :),                                                 &
            real(x_arr(ng : nx_block - ng, :)), real(y_arr(ng : ny_block - ng, :)), real(z_arr(ng : nz_block - ng, :)),            &
            nx_block - 2 * ng + 1, ny_block - 2 * ng + 1, nz_block - 2 * ng + 1, n_blocks_all(ti), n_var, var_dim, format, header)
      endif

!      if (output_method .eq. 3) then
!         call write_vtk_rectilinear (trim(datavis_name) // '.vtk', var_name,                                                       &
!            var(:, ng : nx_total - ng, ng : ny_total - ng, ng : nz_total - ng, :),                                                 &
!            x_base(ng : nx_total - ng), y_base(ng : nx_total - ng), z_base(ng : nx_total - ng),                                    &
!            nx_total - 2 * ng + 1, ny_total - 2 * ng + 1, nz_total - 2 * ng + 1, n_var, var_dim, format, header)
!      endif

!      if (output_method .eq. 4) deallocate(var, x_arr, y_arr, z_arr,                                                               &
!         mesh_name_silo, mesh_name_len_silo, mesh_type_silo, var_name_silo, var_name_len_silo,                                     &
!         level_id_silo, level_data_silo, level_type_silo, level_index_silo, level_data_id_silo)
      if (output_method .eq. 3) deallocate(var, x_arr, y_arr, z_arr)
      if (output_method .eq. 1) deallocate(var, x_arr, y_arr, z_arr)

   enddo

   deallocate (t, dt, n_blocks_all, p_pos,                                                                                         &
               x, y, z, xh, yh, zh,                                                                                                &
               ro, mx, my, mz, bx, by, bz, en, gp, gx, gy, gz,                                                                     &
               roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh,                                                            &
               fgxh, fgyh, fgzh, fggxh, fggyh, fggzh, fmgxh, fmgyh, fmgzh,                                                         &
               fmtxh, fmtyh, fmtzh, jxh, jyh, jzh,                                                                                 &
               var_name, var_dim)

   call mpi_barrier(mpi_comm_world, error)

   if (cpu_id .eq. 0) then
      write(*, *)
      write(*, '(a)') '|--------------------|'
      write(*, '(a)') ' Data manage finished.'
      write(*, '(a)') '|--------------------|'
      write(*, *)
   endif

   call mpi_finalize(error)

!----------------------------------------------------------------------------------------------------------------------------------|
   stop
end program data_manage
